1 Background

Major questions:

  1. Did we induce CTA in red foxes?
  2. How long did it take to induce CTA?
  3. Did we induce CTA in ravens?

2 Data preparation

Firstly, we’re going to install and load our packages. We will use the pacman Package Management Tool which allows us to install and load subsequent packages in a condensed and efficient way.

# Install pacman package
install.packages('pacman')
# Install and load required packages
pacman::p_load(dplyr, ggmap, ggspatial, ggplot2, janitor, lme4, 
               readxl, sf, tidyverse)

Secondly, we’ll read in our data and explore it.

# Read in data
data <- read_excel(path="data.xlsx", sheet="data") %>%
  clean_names()

Let’s view the data.

2.1 Bait take

Now, we’ll wrangle some of our data to get it ready for analyses and plotting.

data <- read_excel(path="data.xlsx", sheet="data") %>%
  clean_names() %>%
  mutate(date = as.Date(date)) %>%
  mutate(bait_take = as.numeric(bait_take)) %>%
  mutate(phase = as.factor(phase)) %>%
  drop_na(bait_take) %>%
  mutate(bait_yes = ifelse(bait_take == 1, "Yes", "No"))

2.2 Weather

The weather data were acquired from the Australian Bureau of Meteorology weather station 70351.

weather <- read_excel(path="data.xlsx", sheet="weather") %>%
  clean_names() %>%
  mutate(date = as.Date(date))

# Combine bait take and weather data into a single df
data <- left_join(data, weather, by=c("date"))

2.3 Species

# Subset to actual instances of bait take by a vertebrate
data <- data %>%
  subset(species != "Lizard") %>%
  mutate(fox_visited = ifelse(fox_visit_1 != 0, "Yes", 
                       ifelse(fox_visit_2 !=0, "Yes",
                       ifelse(fox_visit_3 !=0, "Yes", 
                       ifelse(fox_visit_4 !=0, "Yes", "No"))))) %>%
  
  mutate(fox_visits = ifelse(fox_visit_4 != 0, 4, 
                      ifelse(fox_visit_3 !=0, 3,
                      ifelse(fox_visit_2 !=0, 2, 
                      ifelse(fox_visit_1 !=0, 1, 0))))) %>%

  mutate(raven_visited = ifelse(raven_visit_1 != 0, "Yes", 
                         ifelse(raven_visit_2 !=0, "Yes",
                         ifelse(raven_visit_3 !=0, "Yes", "No")))) %>%
  
  mutate(raven_visits = ifelse(raven_visit_3 !=0, 3, 
                        ifelse(raven_visit_2 !=0, 2,
                        ifelse(raven_visit_1 !=0, 1, 0))))

write.csv(data, "data processed.csv")

3 Analyses

3.1 Bait take

3.1.1 Plots

3.1.1.1 All species

Bait take by treatment

# Read in processed data
data <- read.csv("data processed.csv") %>%
  clean_names() %>%
  mutate(date = as.Date(date))

# Plot by treatment type
bait_treat <- ggplot(data, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("All species") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_treat)

# Save the plot as a jpeg
ggsave("figures/bait_treat.jpeg", height=1500, width=2000, unit="px")

Bait take by treatment and phase

# Read in processed data
data <- read.csv("data processed.csv") %>%
  clean_names() %>%
  mutate(date = as.Date(date))

# Plot by treatment type
bait_treat_phase <- ggplot(data, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("All species") +
  facet_wrap(~phase, scales="free") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_treat_phase)

# Save the plot as a jpeg
ggsave("figures/bait_treat_phase.jpeg", height=1500, width=2000, unit="px")

Percentage bait taken by each species by phase

# Subset to only potential foxes
taken <- data %>%
  subset(bait_take != 0) %>%
  subset(species != "Not taken") %>%
  group_by(species, treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)

# Plot by treatment type
taken_species_phase_hist <- ggplot(data=taken, 
                                        aes(x = species, y = perc, 
                                            col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  facet_wrap(~phase) +
  ggtitle("All species") +
  xlab("Species") + 
  ylab("Percentage of taken baits (%)") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(taken_species_phase_hist)

# Save the plot as a jpeg
ggsave("figures/taken_species_phase_hist.jpeg", height=1500, width=2000, unit="px")

Plot bait take by rainfall

# Plot by rainfall type
bait_rain <- ggplot(data, aes(x = rainfall_mm, 
                                          y = bait_take, 
                                          col=treatment, 
                                          fill=treatment)) +
  geom_smooth() + 
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("All species") +
  xlab("Rainfall") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")
  
# Display the plot
print(bait_rain)

# Save the plot as a jpeg
ggsave("figures/bait_rain.jpeg", height=1500, width=2000, unit="px")

3.1.1.2 Foxes

Percentage fox visits by treatment and phase

# Subset to only potential foxes
fox_stats <- data %>%
  subset(species == "Fox") %>%
  group_by(treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)
  
# Plot by treatment type
visits_fox_histogram <- ggplot(data=fox_stats, 
                              aes(y = perc, x = treatment, 
                                  col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  facet_wrap(~phase) +
  ggtitle("Potential fox visits (including unknowns)") +
  xlab("Species") + 
  ylab("Percentage of fox visits (%)") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(visits_fox_histogram)

# Save the plot as a jpeg
ggsave("figures/visits_fox_histogram.jpeg", height=1500, width=3000, unit="px")
# Subset to only potential foxes
foxes_only <- data %>%
  subset(species != "Raven") %>%
  subset(species != "Unknown")
  
# Plot by treatment type
bait_treat_phase <- ggplot(data=foxes_only, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("Foxes vs no take (excluding unknowns)") +
  facet_wrap(~phase, scales="free") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_treat_phase)

# Save the plot as a jpeg
ggsave("figures/bait_foxes_treat_phase.jpeg", height=1500, width=3000, unit="px")

3.1.1.3 Ravens

Percentage fox visits by treatment and phase

# Subset to only potential foxes
raven_stats <- data %>%
  subset(species == "Raven") %>%
  group_by(treatment, phase) %>%
  summarise(perc = n() / nrow(data) * 100)
  
# Plot by treatment type
visits_raven_hist <- ggplot(data=raven_stats, 
                              aes(y = perc, x = treatment, 
                                  col=treatment, fill=treatment)) +
  geom_col() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  facet_wrap(~phase) +
  ggtitle("Potential raven visits (including unknowns)") +
  xlab("Species") + 
  ylab("Percentage of raven visits (%)") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(visits_raven_hist)

# Save the plot as a jpeg
ggsave("figures/visits_raven_hist.jpeg", height=1500, width=3000, unit="px")
# Subset to only potential foxes
ravens_only <- data %>%
  subset(species != "Fox") %>%
  subset(species != "Unknown")
  
# Plot by treatment type
bait_raven_treat_phase <- ggplot(data=ravens_only, aes(x = as.Date(date), 
                                           y = bait_take, 
                                           col=treatment, 
                                           fill=treatment)) +
  geom_smooth() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust=0.5)) +
  ggtitle("Ravens vs no take (excluding unknowns)") +
  facet_wrap(~phase, scales="free") +
  xlab("Date") + 
  ylab("Bait take") + 
  labs(fill="Treatment", col="Treatment")

# Display the plot
print(bait_raven_treat_phase)

# Save the plot as a jpeg
ggsave("figures/bait_ravens_treat_phase.jpeg", height=1500, width=3000, unit="px")

Bait take by fox visits

# Plot bait take by species visits
bait_visits <- ggplot(data, mapping=aes(x = as.Date(date))) +
  geom_smooth(data, mapping=aes(y = as.numeric(bait_take)), 
              color = "red", fill = "red") +
  geom_smooth(data, mapping=aes(y = as.numeric(fox_visits)), 
              color = "darkorange", fill = "darkorange")  +
  geom_smooth(data, mapping=aes(y = as.numeric(raven_visits)), 
              color = "black", fill = "black")  +
  theme_minimal() + 
  facet_wrap(~treatment) +
  scale_y_continuous(name = "Bait take (%)", limits = c(0, 1), 
                     sec.axis = sec_axis(~., name = "Fox (orange) and raven (black) visits", 
                                         breaks = c(0, 1))) +
  xlab("Date")

# Display the plot
print(bait_visits)

# Save the plot as a jpeg
ggsave("figures/bait_visits.jpeg", height=1500, width=2000, unit="px")

3.1.2 Models

3.1.2.1 All species

Build generalised linear models for our research questions.

# Model for bait take by treatment
summary(glm(bait_take ~ treatment, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment, data = data)
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.66038    0.02223  29.705  < 2e-16 ***
## treatmentTreatment -0.10124    0.03164  -3.199  0.00142 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2357498)
## 
##     Null deviance: 224.02  on 941  degrees of freedom
## Residual deviance: 221.60  on 940  degrees of freedom
## AIC: 1316.1
## 
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and rainfall
summary(glm(bait_take ~ treatment + rainfall_mm, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment + rainfall_mm, data = data)
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.638871   0.023237  27.493  < 2e-16 ***
## treatmentTreatment -0.100640   0.031505  -3.194  0.00145 ** 
## rainfall_mm         0.006263   0.002060   3.040  0.00243 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2337005)
## 
##     Null deviance: 224.02  on 941  degrees of freedom
## Residual deviance: 219.44  on 939  degrees of freedom
## AIC: 1308.9
## 
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and phase
summary(glm(bait_take ~ treatment + phase, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment + phase, data = data)
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.69484    0.05679  12.234  < 2e-16 ***
## treatmentTreatment -0.10146    0.03165  -3.205  0.00139 ** 
## phase              -0.01787    0.02710  -0.659  0.50980    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2358916)
## 
##     Null deviance: 224.02  on 941  degrees of freedom
## Residual deviance: 221.50  on 939  degrees of freedom
## AIC: 1317.7
## 
## Number of Fisher Scoring iterations: 2
# Model for bait take by treatment and phase
summary(glm(bait_take ~ treatment + site, data=data))
## 
## Call:
## glm(formula = bait_take ~ treatment + site, data = data)
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.901039   0.148702   6.059 2.01e-09 ***
## treatmentTreatment -0.102424   0.065102  -1.573 0.116011    
## siteF01             0.046330   0.173125   0.268 0.789062    
## siteF02             0.201385   0.170912   1.178 0.238991    
## siteF03            -0.014532   0.174168  -0.083 0.933523    
## siteF04            -0.157455   0.168687  -0.933 0.350858    
## siteF05            -0.134257   0.170471  -0.788 0.431161    
## siteF06            -0.051039   0.171986  -0.297 0.766718    
## siteF07             0.098961   0.170001   0.582 0.560632    
## siteF08            -0.003736   0.173070  -0.022 0.982781    
## siteF09             0.161306   0.167225   0.965 0.335005    
## siteF10            -0.025888   0.171773  -0.151 0.880239    
## siteF11            -0.003736   0.173070  -0.022 0.982781    
## siteF12            -0.032747   0.169541  -0.193 0.846886    
## siteF13            -0.466256   0.169131  -2.757 0.005956 ** 
## siteF14            -0.560520   0.172711  -3.245 0.001216 ** 
## siteF15            -0.585677   0.170318  -3.439 0.000611 ***
## siteF16            -0.374846   0.168687  -2.222 0.026524 *  
## siteF17            -0.310130   0.170001  -1.824 0.068444 .  
## siteF18            -0.480433   0.171773  -2.797 0.005270 ** 
## siteF19            -0.621452   0.171156  -3.631 0.000298 ***
## siteF20            -0.571342   0.171773  -3.326 0.000916 ***
## siteF21            -0.651039   0.171986  -3.785 0.000164 ***
## siteF22            -0.610447   0.170471  -3.581 0.000361 ***
## siteF23            -0.466256   0.169131  -2.757 0.005956 ** 
## siteF24            -0.578201   0.169541  -3.410 0.000678 ***
## siteF25            -0.453736   0.173070  -2.622 0.008898 ** 
## siteF26            -0.434979   0.171773  -2.532 0.011502 *  
## siteF27            -0.466256   0.169131  -2.757 0.005956 ** 
## siteF28             0.201385   0.170912   1.178 0.238991    
## siteF29             0.196729   0.171156   1.149 0.250693    
## siteF30             0.106147   0.172711   0.615 0.538981    
## siteF31            -0.232064   0.172070  -1.349 0.177788    
## siteF32             0.058162   0.169541   0.343 0.731635    
## siteF34            -0.084329   0.172711  -0.488 0.625480    
## siteF35             0.008052   0.170001   0.047 0.962232    
## siteF40            -0.666907   0.171156  -3.896 0.000105 ***
## siteF48            -0.751374   0.174168  -4.314 1.78e-05 ***
## siteF52             0.104719   0.167760   0.624 0.532646    
## siteF53             0.103414   0.168687   0.613 0.539996    
## siteF54            -0.537745   0.170912  -3.146 0.001708 ** 
## siteF55            -0.091515   0.170949  -0.535 0.592552    
## siteF56             0.008601   0.170471   0.050 0.959774    
## siteF57            -0.526493   0.167171  -3.149 0.001690 ** 
## siteF58            -0.810130   0.170001  -4.765 2.20e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1493372)
## 
##     Null deviance: 223.41  on 937  degrees of freedom
## Residual deviance: 133.36  on 893  degrees of freedom
##   (4 observations deleted due to missingness)
## AIC: 924.16
## 
## Number of Fisher Scoring iterations: 2

3.1.2.2 Foxes

# Do fox visits change with treatment or phase?
summary(glm(fox_visits ~ treatment + phase + 
              rainfall_mm + maximum_temperature_c, data=data))
## 
## Call:
## glm(formula = fox_visits ~ treatment + phase + rainfall_mm + 
##     maximum_temperature_c, data = data)
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)  
## (Intercept)           -0.181054   0.226700  -0.799   0.4247  
## treatmentTreatment     0.022190   0.050005   0.444   0.6573  
## phase                  0.036098   0.053793   0.671   0.5024  
## rainfall_mm            0.007035   0.003418   2.058   0.0398 *
## maximum_temperature_c  0.014803   0.009938   1.490   0.1367  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.5886676)
## 
##     Null deviance: 556.46  on 941  degrees of freedom
## Residual deviance: 551.58  on 937  degrees of freedom
## AIC: 2181.1
## 
## Number of Fisher Scoring iterations: 2

3.1.2.3 Ravens

# Model for visits against treatment and phase
summary(glm(raven_visits ~ treatment + phase + 
              rainfall_mm + maximum_temperature_c, data=data))
## 
## Call:
## glm(formula = raven_visits ~ treatment + phase + rainfall_mm + 
##     maximum_temperature_c, data = data)
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)            0.120731   0.101268   1.192    0.233
## treatmentTreatment     0.016744   0.022338   0.750    0.454
## phase                 -0.002641   0.024030  -0.110    0.913
## rainfall_mm           -0.002171   0.001527  -1.422    0.155
## maximum_temperature_c -0.001225   0.004439  -0.276    0.783
## 
## (Dispersion parameter for gaussian family taken to be 0.1174665)
## 
##     Null deviance: 110.37  on 941  degrees of freedom
## Residual deviance: 110.07  on 937  degrees of freedom
## AIC: 662.88
## 
## Number of Fisher Scoring iterations: 2

3.2 GIS

3.2.1 Map responses

Here we map our responses of bait_take, fox_visits, and raven_visits across Ginninderry Conservation Corridor.

# Use Google API to fetch a base map
# ggmap::register_google(key="[add API key here]", write=TRUE)
map <- get_map(location=c(lon=148.9811, lat=-35.2350),
               zoom=15, source="google", maptype="satellite", crop=FALSE)
# Read in shapefile
gis <- st_read("shapefiles/bait_stations_egg_cta.shp") %>%
  fortify() %>%
  clean_names() %>%
  rename(site=name)
## Reading layer `bait_stations_egg_cta' from data source 
##   `X:\Honours\R2\shapefiles\bait_stations_egg_cta.shp' using driver `ESRI Shapefile'
## Simple feature collection with 43 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 148.9673 ymin: -35.24309 xmax: 148.9927 ymax: -35.22524
## Geodetic CRS:  GDA2020
# Map the bait stations
gin_map <- ggmap(map) + 
  geom_point(data=gis, size=2, aes(x=long, y=lat, col=treatment)) + 
  theme_minimal() +
  xlab("") + ylab("") + labs(col="Treatment")

# Display the plot
print(gin_map)

# Combine the bait take df with the GIS df
data <- left_join(data, gis, by="site")

# Generalised linear model
mod <- glm(bait_take ~ lat + long + elevation, data=data)
summary(mod)
## 
## Call:
## glm(formula = bait_take ~ lat + long + elevation, data = data)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.841e+03  7.170e+02  -6.753 2.56e-11 ***
## lat         -1.071e+01  3.202e+00  -3.345 0.000857 ***
## long         2.997e+01  4.978e+00   6.020 2.51e-09 ***
## elevation   -7.231e-05  6.421e-04  -0.113 0.910369    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1887321)
## 
##     Null deviance: 222.12  on 930  degrees of freedom
## Residual deviance: 174.95  on 927  degrees of freedom
##   (11 observations deleted due to missingness)
## AIC: 1095.7
## 
## Number of Fisher Scoring iterations: 2
??annotation_scale
# Map the bait stations
mod_map <- ggmap(map) + 
  geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
             aes(x=long, y=lat, fill=bait_take)) + 
  theme(axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.ticks.y=element_blank(), 
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(fill="#363842"), 
        strip.text=element_text(colour="white")) +
  facet_wrap(~phase) +
  scale_fill_viridis_c() +
  xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Bait take")

# Display the plot
print(mod_map)

# Map the fox visitation
mod_map <- ggmap(map) + 
  geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
             aes(x=long, y=lat, fill=fox_visits)) + 
  theme(axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.ticks.y=element_blank(), 
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(fill="#363842"), 
        strip.text=element_text(colour="white")) +
  facet_wrap(~phase) +
  scale_fill_viridis_c() +
  xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Fox visits")

# Display the plot
print(mod_map)

# Map the fox visitation
mod_map <- ggmap(map) + 
  geom_point(data=data, shape=21, size=2, stroke=1, col="white", alpha=0.2,
             aes(x=long, y=lat, fill=raven_visits)) + 
  theme(axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.ticks.y=element_blank(), 
        plot.title=element_text(hjust=0.5),
        strip.background=element_rect(fill="#363842"), 
        strip.text=element_text(colour="white")) +
  facet_wrap(~phase) +
  scale_fill_viridis_c() +
  xlab("") + ylab("") + ggtitle("Phase") + labs(fill="Raven visits")

# Display the plot
print(mod_map)

3.2.2 Calculate influence of spatial features on responses

Here we map the bait stations across Ginninderry Conservation Corridor using shapefiles for Ginninderry’s boundaries, paths, roads, and waterbodies.